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Pairing correlations in nuclei play a decisive role in determining nuclear drip-lines, binding ener¬ 
gies, and many collective properties. In this work a new Configuration-Space Monte-Carlo (CSMC) 
method for treating nuclear pairing correlations is developed, implemented, and demonstrated. In 
CSMC the Hamiltonian matrix is stochastically generated in Krylov subspace, resulting in the 
Monte-Carlo version of Lanczos-like diagonalization. The advantages of this approach over other 
techniques are discussed; the absence of the fermionic sign problem, probabilistic interpretation of 
quantum-mechanical amplitudes, and ability to handle truly large-scale problems with defined pre¬ 
cision and error control, are noteworthy merits of CSMC. The features of our CSMC approach are 
shown using models and realistic examples. Special attention is given to difficult limits: situations 
with non-constant pairing strengths, cases with nearly degenerate excited states, limits when pairing 
correlations in finite systems are weak, and problems when the relevant configuration space is large. 
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I. INTRODUCTION 

Pairing correlations are a salient component of the nu¬ 
clear many-body dynamics which has a profound impact 
on most of the nuclear properties, and on the nuclear 
landscape in general. The recently published volume 
Fifty Years of Nuclear BCS [l!] offers a unique overview of 
more than fifty years of research in this area. Advances 
in experimental techniques and emergence of many new 
facilities made it possible to explore nuclear systems at 
the edge of stability. Interest in pairing is reinvigorated 
by extraordinary effects observed recently in near-drip 
line nuclei. This includes the so-called nuclear halo effect 
in neutron rich nuclei, such as the case of borromean nu¬ 
cleus 11 Li which is bound only due to the specifics of the 
pairing dynamics of the two neutrons above the 9 Li core. 
Recent observations of di-neutron decay @1 and other re¬ 
markable manifestations of pairing, seen in structure and 
reactions with exotic nuclei a a, all encourage theoret¬ 
ical effort to be continued. 

The pairing interaction involves pairs of time- 
conjugate single-particle states. We use p\, pk , and hk 
to denote the corresponding pair creation, pair annihila¬ 
tion, and number-of-pairs operators; index k is used to 
identify various distinct pair-states in the system. The 
pairing Hamiltonian of interest is defined as: 

H — 2 ^ ^ Ckfik ^ ' Gkk'PkPk' 7 (f) 

k k,k' 

Here e*, are the single-particle energies and Gkk’ are the 
matrix elements of pairing interaction. In this work we 
limit our discussion to fully paired systems of n pairs in w 
pair-states, which corresponds to 2 n fermions within the 
total particle capacity 2 lo of the valence space. Working 
under the assumption of a fully paired state is completely 
general: any unpaired nucleons remain untouched by the 
Hamiltonian these nucleons effectively block some 
part of the valence space so that the problem is then 
reduced to a fully paired state in a reduced space. Thus, 


the configuration space of interest spans over all ui choose 
n basis states 

|n) = |m,n 2 ,.. .n w ). (2) 

Here we use occupation representation where for each 
pair-state nk = (n|nfc|n) = 1 or 0 depending on whether 
the pair-state is occupied or not. Clearly, the total 
number of pairs n = A n y state can be rep¬ 

resented as a linear combination of the basis states, 

1$) = E n ( n l $ > Im¬ 
pairing correlations have been traditionally explored 
with the help of the BCS theory of superconductivity 
Q. This variational technique, which is formally exact in 
thermodynamic limit, is very well integrated into more 
general mean-held approaches and into techniques be¬ 
yond mean-held. Starting from pioneering works i a, 
the BCS theory has been applied in nuclear physics with 
great success. However, non-conservation of the parti¬ 
cle number and difficulty in handling limits where pair¬ 
ing is weak as compared to the characteristic mean-held 
single-particle level spacing, have proven to be signihcant 
drawbacks in applications of BCS to hnite nuclear sys¬ 
tems [l|; i-El . Over the years a number of remedies have 
been proposed to overcome these drawbacks. For exam¬ 
ple, the issue of the particle number non-conservation has 
been addressed with a variety of techniques proposed in 

Refs. mm. 

With theoretical and computational advances a grow¬ 
ing number of pairing problems in mesoscopic systems, 
such as atomic nuclei, can be treated exactly; thus avoid¬ 
ing the BCS and its drawbacks. There are several ma¬ 
jor groups of exact methods. Symmetry-based algebraic 
methods were introduced by Racah [2CH22I ] even before 
the BCS theory. These methods found wide applicability 
both independently fEi [27| and as components of other 
techniques [El. IE1|. 

Presented more than 40 years ago by Richardson & 
[34j , an exact solution that reduces the pairing eivenvalue 
problem to a set of non-linear equations have been sue- 
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cessfully generalized and applied [HHHI in multiple sit¬ 
uations. Some generalizations and interpretations, such 
as those related to electrostatic analogies [39], are of par¬ 
ticular theoretical interest [ddj . 

Computational advances and iterative sparse matrix 
diagonalization algorithms allowed for direct diagonal- 
ization methods to emerge as extremely simple, stable, 
and robust alternatives |40h42| . Nevertheless, the dimen¬ 
sion of the Hamiltonian matrix in the relevant basis space 
grows exponentially with the number of pairs, eventually 
rendering these methods computationally impractical es¬ 
pecially for model spaces required for problems with con¬ 
tinuum of scattered states. Mote Carlo approach, which 
is the main subject of this work, so far appears to pro¬ 
vide the only reasonable technique that overcomes, in a 
controlled way, the exponentially growing computational 
difficulty. There exist numerous variations of the Monte 
Carlo approach, varying in philosophy and implementa¬ 
tion. Many of these methods can be found in the text¬ 
book [43|. The Shell Model Monte-Carlo [3 [3] and 
Quantum Monte-Carlo involving variational, auxiliary- 
field and Green’s function versions (see review [3]) are 
among well-known successful examples used in low en¬ 
ergy nuclear physics. Random Monte-Carlo sampling, 
either for variational purposes or in order to evaluate 
multi-dimensional integrals, such as those emerging in 
Hubbard-Stratanovich transformation is at the center of 
these techniques. 

The pairing Hamiltonian m is equivalent to a Hamil¬ 
tonian describing w spin-1/2 particles, where one can as¬ 
sume rifc = 1 and 0 for up and down spin orientations, 
respectively, Gamin. This offers opportunities for a 
broad class of Monte-Carlo methods known in spin sys¬ 
tems [13 to be applied. The idea of using the connection 
between spin physics in condensed matter and quasispin 
in pairi ng p roblems was originally explored by Cerf and 
Martin [3 S3- I n their approach the ground state |To) 
is found as an asymptotic state that emerges from an ar¬ 
bitrary initial state $ as a result of evolution along the 
imaginary time: 

|$) ~ e _TBo (T 0 |i , )|4'o}, as r-too. (3) 

In order to propagate the initial state along the imaginary 
time, Cerf and Martin proposed breaking the Hamilto¬ 
nian into the two non-commuting parts Hi and H 2 , cor¬ 
responding to one-body and two-body terms. Then the 
propagation can be done in small steps At using Trotter- 
Suzuki operator decomposition |49l |50|, 

e -Ar(H 1 + H 2 ) = e -Ar^ e -ArH 2 e -Ar^ +Q (4) 

The principal advantage of the technique is that Hi is 
diagonal in the basis states |n) and the corresponding 
exponent can be easily evaluated. Assuming a constant 
pairing, where Gkk' = G, Cerf and Martin proposed to 
evaluate e~ ArH2 stochastically by breaking the exponent 
into a Taylor series and taking advantage of the fact that 
for constant pairing strength the probability of a walk in 


configuration space to have a given number of steps is 
exactly Poissonian. Given that for constant pairing H 2 
can be diagonalized analytically using quasispin algebra, 
it may be possible to use the Trotter-Suzuki propagation 
without involving Monte-Carlo. 

With some degree of success, the method of Cerf and 
Martin was picked up recently by other research groups 
[5ll . I52I ] . The algorithm has so far been applied mainly 
to problems with constant pairing matrix elements. This 
apparent limitation is not well addressed in the literature, 
but appears to be related to unknown quality and relia¬ 
bility of the stochastic evaluation of exponential operator 
of the two-body interaction with non-constant matrix el¬ 
ements. 

In the Configuration-Space Monte-Carlo (CSMC) al¬ 
gorithm presented in this work we completely avoid the 
imaginary time evolution and Trotter-Suzuki decompo¬ 
sition; instead, using a stochastic process of nucleon- 
pair diffusion through the configuration space, we built a 
Krylov subspace which contains the set of lowest eigen¬ 
values. The resulting algorithm can be seen as a Monte- 
Carlo version of the well-known Lanczos algorithm. This 
class of algorithms are often referred to as projector al¬ 
gorithms [43( since repeated application of the Hamilto¬ 
nian operator to a random state eventually amounts to 
the ground state being projected out. Excited states can 
be obtained as well by storing the wave functions and by 
enforcing orthogonality. 


II. CONFIGURATION SPACE MONTE-CARLO 

In this section we present the Configuration Space 
Monte-Carlo method. It should be mentioned that the 
method is generic, and is not limited to pairing, but the 
specifics of the pairing Hamiltonian offer some big advan¬ 
tages, which is discussed in Sec. Ill Cl and demonstrated 
in Sec. IIIII 


A. CSMC formalism 

Let us consider a sequence of states 

|$ l ) = U l |$o) (5) 

which is generated by a repeated application of the 
Hamiltonian H = —V onto a random initial vector |d>o). 
These states span over the Krylov subspace. Eigenval¬ 
ues of the Hamiltonian matrix in this subspace converge, 
after enough iterations, to the eigenvalues (greatest in 
absolute value) of the Hamiltonian in the entire space. 
Since we are interested in the lowest, most negative, 
states it is convenient to carry out this discussion us¬ 
ing V = —H. The repeated application of the opera¬ 
tor V can be written as a summation over all possible 
L + 1 intermediate states which are given by the sets 
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{n}i = {n 0 ,ni,...n L }, 

\*l)= E KM (Ml), (6) 

{"}i 

where the amplitude is 


A ({n} £ ) = (n L |y|n L _i)(n L _i|y|n L _ 2 ) 

... (ni|V|n 0 )(no|$o). (7) 


One advantage of evaluating powers of the Hamiltonian 
operator is that the summation in Eq. © is restricted 
to all possible paths no —> ni • • • —> n £ where each con¬ 
secutive configuration is connected to the previous one 
by the matrix element of the interaction V. Therefore, in 
what follows {n} £ denotes a connected L-step long path 
n 0 -t ni- > n L . 

The path summation can be performed using Monte- 
Carlo sampling. To be more specific, if one generates N 
paths {n}^ = n c ( , s) —> •••—>• n^, labeled here with 

superscript s = 1... N, then 

\*L)*±'t\nP)B({n}U) ( 8 ) 


where 


B{{ n}i s) ) = 


P(Ml } ) 


( 9 ) 


is the amplitude for the s-th random path weighted 
by the inverse of the probability to generate this path 
^({nji )■ I 11 applications where sampling is done with 
uniform probability , P({n}^) = V , the common term 
1/V is just the total number of all possible paths {n} £ 
which equals to the number of terms in the sum in 
Eq. ©. 

Each sampling path can be generated as a random 
walk. The amplitude in Eq. © is subject to the recursion 
relation 


A ({n} L+ i) = (n L+ i|E|n L ) A ({n} £ ), (10) 

where A({n} 0 ) = (n 0 |<f>o). Similarly, the probability 
for a path is the product of probabilities for each step. 
Therefore, starting from the probability to pick the first 
configuration 'P(no) = V ({n}o) , the probability for the 
entire path is generated recursively as 


B (Ml) 


(n £ |E|n £ _i) 
V (n £ _i -)• n L ) 


5 (Ml-i). 


(13) 


The probability distribution for selecting an initial posi¬ 
tion in configuration space 'P(no) and the distribution of 
conditional probabilities V (riL-i —t n £ ) describing the 
direction in which each next random step is to be taken, 
are both arbitrary user-supplied functions. Strategies for 
selecting these functions are discussed in what follows. 

It is important that the probability of taking a certain 
step depends only on the current position and not on 
the preceding history, therefore the process represents a 
Markov chain [43| . The computational implementation of 
the Markov Chain Monte-Carlo methods is a well studied 
subject; see Ref. [43[ and references therein. 

The Configuration Space Monte-Carlo approach, de¬ 
fined by Eq. ©, is implemented using an ensemble of N 
“walkers” starting from configurations no; the initial con¬ 
figurations are generated with the probability distribu¬ 
tion 'P(no). Then each walker independently takes L ran¬ 
dom steps; the probability distribution V (n £ —»• n £+ i) is 
used to generate steps. We envision that each walker car¬ 
ries a “bag” B that is initialized and modified along the 
path following Eqs. m and m■ Contributions from 
the bags of all walkers arriving to a given configuration 
n £ comprise the component (n £ |$ £ ) as shown by Eq. 
©. 

As the most straightforward application of the method, 
one could assume the probabilities for steps in all “di¬ 
rections” to be equal, then the conditional probabil¬ 
ity V (n —» n') depends only on the initial configuration 
n and the inverse of it equals to the number of con¬ 
figurations connected to n. In most cases the number 
of connected configurations is the same for all states, 
which makes the conditional probability for each step 
being an absolute constant, i.e., independent of initial 
and final positions. For example, for any paired con¬ 
figuration with n pairs and w pair-spaces the pairing 
Hamiltonian can generally move one of the n pairs onto 
one of the w — n + 1 unoccupied pair-states (this in¬ 
cludes diagonal move back to the same pair-state). Thus, 
in the pairing case, for equiprobable steps the condi¬ 
tional probability becomes a configuration-independent 
constant V (n —>• n') = (n(w — n+1)) -1 and the result¬ 
ing random paths are all generated with equal probabil¬ 
ity. This amounts to uniform Monte-Carlo sampling of 
terms in sum ©. 


B. Importance sampling 


B (Ml+i) = B(n L ->■ n L+ i) V ({n} £ ); (11) 


here V(n Q —>• n £+ i) is the conditional probability to 
move to configuration n £+ i given the current position 
at n £ . Thus, while going along a random path the coef¬ 
ficients B are generated recursively, 


B{ Mo) 


(no|$p) 

Vino) 


and 


( 12 ) 


Uniform sampling is convenient and effective when con¬ 
tributions from most paths are nearly equal; constant- 
strength pairing Hamiltonian discussed by Cerf and Mar¬ 
tin in Refs. [13, S3 is a good example of this situa¬ 
tion. However, sampling uniformly can be extremely in¬ 
effective if certain amplitudes A (|n} £ ) are very small or 
equal to zero; importance sampling can be introduced as 
a remedy. In the CSMC the contributions from different 
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sampling paths can be made comparable in magnitude 
if steps are generated with probabilities proportional to 
the magnitude of the corresponding matrix elements, 

V (riL-i -A n L ) ex |{nL|V r |n i _i}|. (14) 

This way the scaling factor in Eq. (1131) would not depend 
on the direction of the step. It should be emphasized, 
that satisfying the proportionality m exactly, which 
may be computationally expensive, is not necessary. Any 
probability distribution that in some general way follows 
the distribution of the matrix elements is sufficient. 

The approach described here, referred to as Configura¬ 
tion Space Monte Carlo, allows one to build stochastically 
the Krylov subspace and find the eigenstates and eigen¬ 
values of the Hamiltonian using steps similar to those in 
Lanczos approach. Clearly, the method is applicable to 
any Hamiltonian; however, different signs of matrix ele¬ 
ments (n£_j_i|K|ni) and thus different signs of the am¬ 
plitudes can lead to poorly convergent sums. This issue, 
commonly known as the Monte-Carlo sign problem, is not 
present in applications of the CSMC to pairing problems 
that are discussed next. 


C. Features of the pairing Hamiltonian 

Let us summarize some of the important features of 
the pairing problem that boost the effectiveness of the 
CSMC method. 

(i) For fully paired systems the diagonal pairing matrix 
elements Gkk are equivalent to the single-particle ener¬ 
gies. Thus, by redefining the diagonal pairing matrix 
elements as Gkk -A Gkk — e/c/2, the pairing Hamiltonian 
can be written in the following form 

V = —H, where V = Gkk'PkPk' ■ (15) 

k,k' 

(ii) The pairing interaction is attractive. Therefore, with 
the proper choices of phases all off-diagonal matrix ele¬ 
ments Gkk’ can be made non-negative. Without any loss 
of generality the single-particle energies can be measured 
relative to some chemical potential /x; and the constant 
/x can be selected so that all diagonal many-body matrix 
elements (n|V|n) are also positive. 

(iii) The nucleon pairs are the only degrees of freedom, 
and the entire dynamics is represented by the “hopping” 
of the pairs between available states. Each pair hopping 
leads to a step in configuration space where 

(n'|V|n) = (...n k = 1, ...n k '= 0, ...|V|. ..n/ c =0, ...n k >= 1, ■■■) 

= G kk ' > 0 if k ± k', 

and (n|U|n) = G^n*, > 0 for the diagonal. For any 

initial configuration there are n(u> — n + 1) different final 
configurations that can be reached in one step. 


(iv) Given that the matrix elements of V are all positive, 
the Monte-Carlo sign problem does not appear. 

(v) Positive matrix elements of V imply that if in the ini¬ 
tial wave function $o all components (n|$o) > 0: which 
can always be accomplished by defining phases of the 
basis states |n), then all components of any are non¬ 
negative, that is, (n|$i) > 0 for any L and any n. This 
also allows one to introduce a linear L\ norm 

( 16 ) 

n 

(vi) Asymptotically, as L — > oo, 

|*l> - (-£o) i ('Lo|$o}|*o>, (17) 

where |T 0 ) is the ground state wave function and E 0 is 
the ground state energy. Since Eg < 0, in fact for the 
negative-definite Hamiltonian, all eigenvalues are nega¬ 
tive, and the phase of the ground state wave function 
can be selected so that (H> 0 |<H 0 ) > 0, and all components 
of the ground state wave function are also non-negative: 
(n|T 0 ) > 0. 

(vii) Given that all many-body states that span the 
Krylov subspace have positive-definite amplitudes rel¬ 
ative to the basis states |n), these amplitudes can be 
treated as probabilities. Therefore in the ideal limit of 
the importance sampling Monte-Carlo, when Eq. m is 
satisfied exactly, all walkers’ bags are equal. In this limit 
the number of walkers arriving to a certain many-body 
configuration n on step L is proportional to (n|4>i). Lin¬ 
ear norm can be used to normalize the wave function. 


III. NUCLEAR PAIRING WITH 
CONFIGURATION SPACE MONTE CARLO 

In what follows we demonstrate the CSMC applica¬ 
tions to pairing problems. We organize our presentation 
by progressing from simple to more elaborate applica¬ 
tions, we use model and realistic examples to highlight 
the CSMC and to address technical details. 

In the following subsections A-D we consider a system 
consisting of ui double-degenerate equally-spaced single 
particle orbitals, = e k with k = 0,1... uj — 1, where 
the single-particle level spacing e defines the unit of en¬ 
ergy. This system, often referred to as a picket-fence 
or ladder model, is commonly used for testing various 
approaches to pairing [42| . The model has a minimal 
symmetry, time-reversal only, making it the most com¬ 
putationally challenging one. For the set of studies using 
the ladder model we assume constant pairing interaction, 
where Gkk' = G in Eq. ©. This choice is not related to 
any limitation, this merely minimizes the number of pa¬ 
rameters and is convenient for comparison with numerous 
previous studies where the same model was used. For the 
half-occupied ladder model the critical BCS strength is 
approximately G cr = e/ln(2w), see Ref. [Hi]. 


5 


Realistic examples with non-constant pairing strength 
and large scale application are discussed in subsections 
E and F. 


A. Linear norm 


We start with a very simple and quick technique that 
can be used to determine the ground state energy and 
requires no storage for wave functions. Despite certain 
limitations, the method is elegant, simple in implementa¬ 
tion, very computationally efficient, and is an important 
component in the general approach. 

In the implementation of the CSMC through multi¬ 
ple walkers in configuration space the construction of the 
wave function |d>i) is the most challenging task because, 
according to Eq. (j5]| , it requires organizing walkers based 
on their arrival locations. Even if performed in parallel, 
this is still a daunting task when the number of contribut¬ 
ing configurations becomes large. As we show next, for 
certain observables, and for the ground state energy in 
particular, this task can be avoided; thanks to the prop¬ 
erties of pairing interaction outlined in Sec. Ill Cl 

According to Eq. (j8]) the average of all bags gives the 
linear norm (referred to as C\) of the wave function 


!|>({n}«) « = IM; ( 18 ) 


computing the bag average is a simple and fast operation. 

Following Eq. (ED. the bag averages for two consecu¬ 
tive values of L as L —> oo give an estimate for the ground 
state energy as 


Eq — Eo(L) = 


En( n l^+l) 

E„(n|$L> 


ll<M ' 


(19) 


Clearly, any procedure on the Hamiltonian leading to the 
ground state can be subjected to the linear norm. For 
example, using projection ([3]) one could evaluate energy 
in the limit r-> oo as 


E 0 ~ E 0 (t) 


\\He- rH * 0 \\ 

||e-^$ 0 || 


( 20 ) 


where exponents are evaluated in CSMC approach using 
Taylor series 


OO £ 

\\e~ TH M = Elfll^ll- ( 2 D 

L =0 Ll ' 

Obviously, it is possible to compute the linear norm for 
any operator ||0$|| = JI n (n|0|$), but unfortunately in 
most situations this linear norm does not have a trans¬ 
parent physical meaning. 

This quick linear-norm-based technique for evaluating 
the ground state energy within CSMC algorithm is il¬ 
lustrated in Fig. [TJ and is compared to the general ap¬ 
proach discussed in the following subsection. A ladder 


model with ui = 18 levels and n = 9 nucleon pairs, where 
G = 1, is used in this example. Two different starting 
wave functions |<I>o) are considered. In the Fermi state all 
lowest single-particle levels are occupied up to the Fermi 
surface, 


l^o Fermi) ) = II p!|°> = I 1 • • ■ !. 0. 0 • ■ ■ )• (22) 

k=l " v ' 

n spaces 

The Fermi state is an exact ground state of non¬ 
interacting fermions (G = 0 limit). 

The BCS solution offers a second convenient starting 
wave function n£=i(«* + Ufcp|.)|0). In our applications 
it is projected onto an appropriate number of particles; 
therefore is defined via amplitudes as 

UJ 

(n|$o BCS) ) = II ( UkSn ^o + Ufc<W)- (23) 

fc= l 

The coefficients Uk and vk are determined by solving the 
usual BCS equations. Our procedure does not require 
the starting wave function to be normalized. Given a 
product form of Eq. (l23l) . it is efficient to generate the 
BCS based initial state stochastically by selecting |no) in 
a process where each randomly selected state k is chosen 
to be occupied or empty with a probability proportional 
to the corresponding Vk and Uk', the process is stopped 
once a desired number of occupied states given by the 
total number of particles is reached. It is important that 
variations in implementation of the projected BCS or not 
following Eq. (l23l) exactly, are not essential since the 
starting wave function can be arbitrary. 

In Fig. [T]( a) the convergence of energy as a function of 
L, following Eq. cni), is shown for the two initial wave 
functions. The method just outlined is based on eval¬ 
uation of the linear norm, the sum of all walkers’ bags, 
and since the actual wave functions are never constructed 
these are labeled as “no-wf” in Fig. [T] In panel (b) the 
convergence is shown as a function of the imaginary time 
r using Eq. (l20l) . The common energy scale is used in 
both panels and the energy obtained from BCS approach 
and from the exact diagonalization of the pairing Hamil¬ 
tonian are shown with horizontal grid lines. The magni¬ 
fied energy scale used here allows one to clearly see the 
difference between BCS and the exact solution. 

As appropriate in a variational technique, the BCS 
ground state energy is above the exact one. However, 
the estimates for the ground state energy, using the lin¬ 
ear norm, approach the exact value from below. This 
feature, as discussed in Sec. IIIIBI is used for providing a 
lower bound for the ground state energy estimate. 

In order to compare the projection with power func¬ 
tion in Eq. USD, and using the exponential m Eq. (jldlj) , 
Fig. [T](b) includes an additional L scale shown at the top. 
The quantity L is defined as the average number of steps 
that needs to be taken by walkers in order for the se¬ 
ries (1211) to converge for a given imaginary time r. While 








BCS 

Exact 



FIG. 1: (Color online) The half-occupied ladder model with cj = 18, n = 9 and G = 1 is used to show the convergence of ground 
state energy using an approach based on the linear norm. Left panel (a) shows the method based on ground state projection 
with power law Eq. urn and right panel (b) shows projection using imaginary time evolution Eq. (12011 . In both panels BCS 
and exact values of energy are shown with horizontal grid lines; the four curves represent two initial states and two methods: 
with and without wave functions being built. 


both panels (a) and (b) look similar, using the exponen¬ 
tial as a projector is more computationally expensive as 
it requires almost three times as many steps. 

The use of exponent to project a ground state does not 
provide any additional numerical stability; fluctuations 
at remote times, in cases with no wave function, are seen 
in both panels of Fig. [T] These fluctuations are removed 
by reconstructing wave functions at certain steps; the 
corresponding curves in Fig. |T] are labeled with “wf”. 
The origin of these fluctuations and error analysis are 
addressed next. Since the exponential projection using 
imaginary time is deemed to be less effective we will not 
discuss it any further. 


B. Error and convergence control 

In the CSMC there are generally two kinds of errors. 
The first one is the statistical error that emerges as a 
result of stochastic evaluation, for example, estimating 
wave functions using Eq. ([5]) or evaluation of the linear 
norm in Eq. m- The second error is associated with the 
algorithm used to obtain physical quantities of interest; 
for example, in projection technique this concerns the 
quality of approximation Eq ( L ) ss Eq in Eq. (fl9l) . In this 
subsection we examine both of these errors and methods 
of their control. 

The Central Limit Theorem (CLT) is at the core of 
statistical error control. It is usually expected that as 
the number of samples, N, grows the associated stan¬ 
dard deviation a of the ensemble average goes down as 
a oc 1/y/N. However, in CSMC the main disadvantage 


of independent walks is that the variance grows exponen¬ 
tially as the path length increases. Therefore, for large 
number of steps the average of bags in Eq. ca is hard 
to evaluate because the distribution of bags becomes too 
broad. This is the cause of fluctuations seen in Fig. m 
at large L. 

Let us analyze this problem. Consider an ensemble of 
all L-step bags for all possible paths {Bl}, let ct 2 {Bl} 
be its variance and Bl its mean. According to Eq. (fl8l) 
Bl = ||$l||- As proved earlier, all bags are positive mak¬ 
ing the coefficient of variation Cv{Bl} = ct{Bl}/Bl an 
appropriate measure of relative error. Indeed CLT im¬ 
plies that with N estimates of energy using Eq. (flT)l) the 
relative error is 

AE 0 (L)/E 0 (L) » J=C v{B l+1 }. (24) 

The problem with divergent behavior of Cv{Bl} as a 

(s) 

function of L arises due to B L for each walker s being a 
product of matrix elements weighted by the correspond¬ 
ing probability, see Eq. m- The product of a large 
number of random matrix elements is poorly behaved. 
Let us assume that c gives the coefficient of variation for 
all possible matrix elements weighted by chosen probabil¬ 
ities, then each term in the product (fl3l) has a coefficient 
of variation, 


Cv 


(n L \V\n L -i) \ = c2 
V (n L _i -*• n L ) J ~ 


(25) 


Then the standard statistical treatment of a product 
leads to 


Cv{B l } = \J {l + c 2 ) L — 1; 


(26) 
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this simple form is obtained under the assumption of uni¬ 
form initial distribution in Eq. m- With the exemption 
of some special cases, c > 0; in most situations of inter¬ 
est c ss 0.5, therefore C v{Bl} grows exponentially with 
L. Thus, it is practically impossible to compete with the 
exponentially increasing variance of the distribution by 
increasing the number of walkers. 

In Fig. @ the behavior of C v{Bl} as a function of L 
is shown for the half-occupied 18-level ladder model. The 
results for two different starting states are shown. The 
exponential divergence in Eq. (1261) pertains to the sit¬ 
uation involving independent walkers, where the actual 
wave function is not obtained. The two corresponding 
curves, labeled with “no-wf”, both display the same ex¬ 
ponential divergence with c « 0.42 which corresponds to 
asymptotic behavior, Cv{Bl} oc 1.086 l . 



L 


FIG. 2: (Color online) The half-occupied ladder model with 
n = 9, uj = 18, and G = 1. Coefficient of variation Cv{Bj,} is 
shown as a function of L. The figure includes four curves with 
two possible initial states: (a) Fermi state 3>( Ferml ) Eq. (1221) 
and (b) constant-component state where (n|<f>Q Const ^) = 1; and 
for calculations with and without wave functions being ob¬ 
tained. 

The limitation on the number of independent steps is 
relatively easy to overcome. The exponential growth of 
Cv{Bl} is usually weak, and in most cases, such as the 
example in Fig. [I] no problems emerge for L less than 30 
or 50. Moreover, the choice of probabilities that follows 
importance sampling in Eq. m would lead to c = 0. 
Practically, numerical noise never allows one to reach 
this ideal limit but the the growth of variance can be 
delayed. In addition to that, with a good initial wave 
function, such as the one from BCS theory, the conver¬ 
gence is reached in a few steps, before the onset of sta¬ 
tistical problems; see example in Fig. [T] 

Preventing walkers from taking long independent 
walks, by combining them in wave functions after a cer¬ 
tain number of steps with Eq. ©, allows one to avoid the 


problem completely. This is demonstrated in Fig.[T]with 
curves labeled “wf”. The intermediate summation at a 
moment when bags are combined, due to the CLT, pre¬ 
vents an exponential increase of the variance. In the im¬ 
plementation of the CSMC the statistical error is tracked 
by controlling the coefficient of variation in the bags as 
L is increased. This allows one to apply a computation¬ 
ally expensive procedure of reconstructing the wave func¬ 
tion only when necessary, typically once in every 5-20 
steps. At a moment when the full wave function is built, 
the convergence of the projection technique (the second 
kind of error) can be assessed using the usual square, £2, 
norm. 

The second kind of errors, which is convergence 
Eq(L) —> Eq in these examples, is a part of any itera¬ 
tive diagonalization technique, such as Lanczos or David¬ 
son algorithms; and it has been well studied in the past. 
However, the specifics of the pairing problem described 
in Sec. Ill Cl and the use of the linear norm allow one to 
place exact upper and lower limits on the value of energy. 

The convergence of the projection algorithm is exam¬ 
ined in Fig. [3] Here, using the same half-occupied ladder 
model with uj = 18, we show deviation of the predicted 
energy from the exact value as a function of the number of 
steps. Three curves, that are essentially indistinct, show 
Eq(L) — Eq where Eq(L) was evaluated using the linear 
norm £1, Eq. (fl^l) . The three sets of results are obtained 
by evaluating |$^) exactly with matrix-vector multiplica¬ 
tion (dotted line); using CSMC with wave function being 
reconstructed at each step (dashed line); and without the 
wave function using bag average in Eq. (fT51) (solid line). 
The slight difference between exact and CSMC results 
is only due to an intermediate shift by chemical poten¬ 
tial. These three curves approach the exact energy from 
below, which is a distinct property of the C\ norm. 

The other two, nearly indistinct, curves show the con¬ 
vergence of energy evaluated using the traditional square 
norm, labeled as £2 


Eo(L) 




(27) 


The £2 norm can be used only when the wave function is 
available, for that reason only the curves for exact (dash- 
dot) and CSMC with wave function (short dash) appear 
in Fig. [31 Naturally, the expectation value of the Hamil¬ 
tonian in Eq. (1271) is subject to the variational principle, 
and all curves with £2 norm approach ground state en¬ 
ergy from above. Thus, the estimates using £1, Eq. m, 
and £2, Eq. (1271) norms give the lower and upper bounds 
for the value of ground state energy. 

To summarize, in our algorithm we rely on computa¬ 
tionally inexpensive independent propagation of walkers 
in configuration space until the coefficient of variation 
of their bags exceeds some critical value. At that mo¬ 
ment the full wave function is reconstructed and is used 
to evaluate energy from Eq. (EH) and all other operators 
of interest. The combination of energy estimates from 
linear and square norms give lower and upper bounds for 












L 

FIG. 3: (Color online) The half-occupied ladder model with 
n = 9, w = 18, and G = 1. Deviation of the energy es¬ 
timates using linear C\ and square £2 norms from Eqs. m 
and (E3. respectively, is shown as a function of L. The curves 
correspond to exact, CSMC with and without wave function 
reconstruction. 

the actual value of energy. If the desired convergence is 
not reached, the process is continued starting from the 
current wave function. 


C. Weak pairing limit 

As mentioned earlier, superconducting paired states in 
small systems face a lot of competition from other in¬ 
coherent interactions as well as from the single particle 
shell structure. Thus, relatively weak and fragile super¬ 
conducting states is one of the distinct characteristics of 
pairing in nuclei. Unfortunately, the BCS theory is not 
designed to work in this limit, and having the CSMC as 
a computationally inexpensive alternative is one of the 
main motivations of this work. In Fig. [4] we demonstrate 
the effectiveness of CSMC in the limit of weak pairing 
using our half-filled 18-level ladder model. In the limit 
when the pairing strength G = 0, the system settles in 
the Fermi state with 9 lowest double-degenerate single¬ 
particle states being occupied. As soon as G > 0, pair 
excitation promotes particles up, and the occupation of 
the upper 9 levels becomes non-zero. For very strong 
pairing, G>f, the limit of degenerate model with equal 
occupancy of all states is reached. This limit leads to 
half of the 18 particles being on lower 9 levels and half 
on the upper ones. 

In Fig. |4] the net occupation of the upper 9 levels as a 
function of G is shown. This plot includes results from 
BCS, CSMC (labeled as “MC”) and exact diagonaliza- 
tion. While on a large scale all results are similar, in 


the region of low pairing strength, which is shown in in¬ 
set, the well known problem with BCS solution, shown 
in dotted (black) line, is noticeable. At the same time, 
exact and CSMC results are indistinct; dashed (crimson 
color) line goes right on top of the solid (sea-green color) 
line. This test illustrates that the CSMC is well suited 



G 


FIG. 4: (Color online) The half-occupied ladder model with 
uj = 18 and G = 1. The net occupancy (total number of 
particles) on the upper 9 orbitals is shown as a function of 
the pairing strength. The weak pairing limit is magnified in 
inset. The three curves correspond to BCS solution, CSMC 
(MC) solution, and the exact solution by means of diagonal- 
ization. The CSMC and exact results are indistinct and the 
corresponding curves are overlaid. For the CSMC solution we 
used N = 7.5 x 10 5 walkers, limiting the number of indepen¬ 
dent steps to five. 

for all limits of pairing strength. Moreover, in the limit 
of weak pairing, the computational effort in CSMC is 
reduced as the contribution from rare excursions above 
Fermi surface can be easily evaluated with importance 
sampling. 

D. Excited states 

Obtaining excited states with CSMC is more compu¬ 
tationally difficult. One can no longer use a linear norm 
since all amplitudes cannot be positive definite simulta¬ 
neously; that is, statement (vi) in Sec lII Cl is not valid 
for excited states. Therefore, the usual quadratic £2 
norm has to be used and the bag values can be neg¬ 
ative. Nevertheless, features (i)-(v) in Sec lII Cl remain 
valid and useful. In particular, since the matrix elements 
of V are positive definite, the importance sampling is still 
an effective strategy and the signs of bags are not altered 
by repeated application of the Hamiltonian which cur¬ 
tails the typical MC sign problem. Similar, to Lanczos 
technique, the CSMC approach requires orthogonaliza- 
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tion, therefore the wave functions have to be built each 
time the orthogonalization is to be performed. The need 
for orthogonalization limits the number of independent 
steps, which is the main reason for higher computational 
demand. 

In Fig. [5] we show the CSMC applied to the study of 
excited states in the same half-occupied 18-level ladder 
model. The ladder model example is particularly chal¬ 
lenging since the density of states above the gap is high. 
In this model the level spacing between the ground and 
first excited state, which is about twice the BCS gap, is 
E\ — Eg ft 14.6 (in units of level spacing e = 1). At 
the same time, the spacing between the following states 
E 2 —E 1 ft 0.3 is very small. Moreover, the second excited 
state is double-degenerate. 



L 

FIG. 5: (Color online) The half-occupied ladder model with 
co = 18 and G = 1. Convergence of CSMC to ground state, 
to first excited state and to double-degenerate second excited 
state is shown. 


E. Pairing in Sn isotopes 

In order to illustrate the CSMC algorithm in a realistic 
case where pairing matrix elements are not all equal to 
a constant, we consider isotopes of tin. The role of pair¬ 
ing in 100 ~ 132 S n isotopes has been extensively explored 
in the literature [l0|, HH US [Hi . Apart from questions 
of scientific interest such as pairing matrix elements and 
their connection to superconducting state in infinite mat¬ 
ter, near constancy of the excitation energy of the lowest 
2 + states, and unexplained behavior of electric quadru¬ 
ple transition rates, the tin case emerged as a benchmark 
for computational techniques. In Tab. U we present com¬ 
parison of energies and occupation numbers for 116 Sn, 
118 Sn, and 120 Sn. The model space here includes five 
single-particle levels (with total co = 16), their energies 


and spins are listed in the first two columns of Tab. Q] 
the matrix elements are taken from from the G-matrix 
calculation in Ref. [54|, the values can be found in Ta¬ 
ble 1 of Ref. [T3|. The results in Tab. Q] show expected 
level of agreement. With increased computational effort, 
mainly using large number of walkers, any desired level 
of precision can be obtained; our goal here was to use 
minimal effort and to solve the pairing problem with a 
precision that exceeds any practical need, which is set to 
be 5 keV uncertainty for energy and 0.01 for occupation 
numbers. 


F. Large scale model 

As a final illustration of the CSMC algorithm we ex¬ 
plore a model of the 24 O nucleus intended to reflect the 
nature of pairing correlations in a system containing both 
bound states and a continuum of scattering states. Our 
main goal is to demonstrate the capabilities of our al¬ 
gorithm while addressing the problem of pairing in con¬ 
tinuum qualitatively. Quantitative studies require good 
knowledge of the effective interaction Hamiltonian; con¬ 
struction of this Hamiltonian is outside the scope of this 
presentation. 

For our study we select the Woods-Saxon potential 
with parameters from Ref. [55| to model the mean field of 
weakly-bound 24 0 nucleus. We discretize this potential 
using a large quantization-box of size 500 fm. This al¬ 
lows us to generate a dense continuum of states. We limit 
scattering states by about 8 MeV of energy, which leads 
to co of about 100. For the pairing interaction between 
neutrons we use a density dependent contact interaction 
from Refs. [56i fbijf 

V(r, r') = — G 0 ^1 — <5(r — r'). (28) 

Here p(r)/po is the nucleonic density expressed relative 
to the saturation density. This quantity is assumed to be 
given by the Woods-Saxon form factor. The density de¬ 
pendence of pairing is controlled by a parameter p which 
is selected as p = 0.5. Following Ref. ;56i, H}, we also 
introduce a momentum cut-off function, that gradually 
reduces the pairing matrix elements to zero for scattering 
states at energies above 5 MeV, the diffuseness parameter 
of the cut-off function is 0.5 MeV, see also Ref. [59j . 

For our example we assume an inert 16 O core which 
leaves two bound S 1/2 and d 5 j 2 valence single-particle 
states. Therefore, the bound states can accommodate 
n = 4 pairs of valence neutrons in 24 O. The pairing 
matrix elements involving these states are known from 
the phenomenological shell model Hamiltonian in Ref. 
f60| . Following previous studies, we adopt the value of 
Go = 1 GeV-fm 3 for treating pairing interaction involv¬ 
ing the continuum of scattered states. This value is also 
consistent with t he p airing strength in phenomenologi¬ 
cal Hamiltonians [60j. Here we limit our consideration 
to s-wave single-particle continuum. Due to centrifugal 
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llb Sn 

Exact 

CSMC 

ii8 Sn 

Exact 

CSMC 

i2U Sn 

Exact 

CSMC 

j 

.Eo(MeV) 

-153.766 

-153.765 

-170.115 

-170.113 

-185.945 

-185.945 

5/2 

-9.736 

4.99 

4.99 

5.13 

5.13 

5.25 

5.25 

7/2 

-8.957 

5.88 

5.89 

6.28 

6.27 

6.60 

6.60 

1/2 

-7.302 

0.67 

0.67 

0.76 

0.76 

0.86 

0.86 

3/2 

-7.634 

1.18 

1.17 

1.63 

1.63 

2.08 

2.09 

11/2 

-7.544 

3.29 

3.28 

4.21 

4.20 

5.21 

5.21 


TABLE I: Comparison of exact and CSMC results for selected isotopes of tin. After header, first row shows comparison of 
energies, the remaining five rows show occupation numbers for five singe particle states. The calculation of energies is done 
with N = 5 x 10 6 walkers using linear norm. The final error is about 5 keV. 


barrier the overlap between bound and unbound d-wave 
states is small; this inhibits virtual pair excitations to 
d-wave states in the continuum. 

Our goal in this investigation is to estimate the role 
of continuum. We do this by comparing full calculation 
with the one where the continuum is ignored. In Fig. [G] 
the change in the ground state energy AEq is shown as a 
function of the single-particle energy e of the Si/ 2 state. 
We present two different cases. In case (a) the parame¬ 
ters of the pairing Hamiltonian, which includes the single 
particle energies and pairing matrix elements, are first 
evaluated with a realistic choice of Woods-Saxon param¬ 
eterization for 24 O and then the e which corresponds to 
S \/2 state, is varied while all other parameters remain un¬ 
changed. In the self-consistent case (b) the depth of the 
Woods-Saxon potential is varied which moves the Si/ 2 
state, and each time a new configuration space Hamilto¬ 
nian matrix is calculated and studied. 

Let us summarize this study. First, the correction from 
pair excitations into continuum appears to be relatively 
small, here it is of the order of one kilovolt, for all rea¬ 
sonable choices of pairing strength Go the effect is not 
expected to exceed a few tens of kilovolts, see Ref. 5Sf. 
The smallness of the effect does not seem to contradict 
observations. So far there has been no significant near¬ 
threshold discontinuity observed in nuclear structure that 
can be attributed to two-body decay or to pair excita¬ 
tions. The decay of 26 0 is observed to be very slow, see 
discussion (6l|, which through dispersion relations indi¬ 
cates weakness of the continuum coupling. 

Second, as expected, the effect increases sharply as the 
bound state approaches the continuum threshold. This 
is similar to the results known for single particle states, 
while the exact near threshold behavior is defined by the 
phase space volume, see Refs. [H, ! 63i]. 

Third, the difference between the two models high¬ 
lights the importance of halo phenomenon and its proper 
treatment. In model (b) the wave function of the 
single-particle S\/ 2 state spatially extends as its energy 
approaches the threshold. This facilitates pairing in 
the continuum and the resulting effect is significantly 
stronger than that in model (a) where the spatial struc¬ 
ture of the single particle wave function was not modified. 

Finally, we find that this example successfully demon¬ 


strates the power of the CSMC method. 



FIG. 6: (Color online) The energy correction due to the in¬ 
clusion of continuum states for a pairing model with lo = 100 
and n = 4. The correction, AEo, is the difference between 
the CSMC result with continuum states and an exact answer 
for a model including only the two bound states. The e along 
the lower axis is the energy of the s-wave bound state as it 
is moved closer to the continuum. The Monte Carlo error is 
negligible. 


IV. SUMMARY 

In this work we put forward a new Configuration Space 
Monte Carlo method for solving the many-body pairing 
problem in finite systems. Unlike previous Monte Carlo 
techniques that deal with pairing interaction in a way 
similar to MC methods in physics of spin systems, our 
approach does not use evolution in imaginary time, does 
not need Trotter-Suzuki propagator breakup, and does 
not depend on the pairing matrix elements being con¬ 
stant. We propose to evaluate Hamiltonian and other 
observables by stochastically evaluating the correspond¬ 
ing operators in the Krylov subspace spanned by states 
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formed as a result of powers of Hamiltonian acting on an 
arbitrary initial state. States in the Krylov subspace are 
evaluated using random walks in the many-body config¬ 
uration space. Importance sampling is used to effectively 
probe components of the wave functions. We emphasize 
several important features of the pairing Hamiltonian, 
that make the MC approach appealing. In particular, 
we stress boson-like behavior of nucleon pairs, absence 
of the fermion sign problem, potential for probabilistic 
interpretation of transitions in configuration space, and 
probabilistic interpretation of ground state amplitudes. 
In addition to traditional quadratic quantum mechanical 
norm, probabilistic interpretation allows us to use a lin¬ 
ear norm. We demonstrate that the approach based on 
the linear norm is computationally efficient, is perfect for 
parallelization, and provides effective methods for control 
of errors of both stochastic and non-stochastic origins. 

The workings of the CSMC method are demonstrated 
with several examples. With a classic ladder model we 
demonstrate convergence using several variations of the 
method; we discuss errors and present effective means of 
their control. The effectiveness of CSMC in small sys¬ 


tems where pairing can be effectively weak is shown; the 
CSMC in its most complete form is used for obtaining 
degenerate and nearly-degenerate excited states in the 
ladder model. 

As a realistic example, we use isotopes of tin which 
represents another well studied classic case of pairing 
in nuclei. The energies and occupation numbers in this 
non-constant pairing example are consistent with exact 
results. Large-scale study of pairing correlations is illus¬ 
trated using a model of 24 0 that includes a continuum of 
scattering states. While our last example is still far from 
realistic, it highlights the effectiveness of CSMC, and sug¬ 
gests an arena for future applications of the method. 
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